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Abstract. The fractal dimensions and the percolation exponents of the geometrical 
spin clusters of like sign at criticality, are obtained numerically for an Ising model with 
temperature-dependent annealed bond dilution, also known as the thermalized bond 
Ising model (TBIM), in two dimensions. For this purpose, a modified Wolff single- 
cluster Monte Carlo simulation is used to generate equilibrium spin configurations on 
square lattices in the critical region. A tie-breaking rule is employed to identify non- 
intersecting spin cluster boundaries along the edges of the dual lattice. The values 
obtained for the fractal dimensions of the spanning geometrical clusters Dc, and their 
interfaces D/, are in perfect agreement with those reported for the standard two- 
dimensional ferromagnetic Ising model. Furthermore, the variance of the winding 
angles, results in a diffusivity k = 3 for the two-dimensional thermalized bond Ising 
model, thus placing it in the universality class of the regular Ising model. A finite-size 
scaling analysis of the largest geometrical clusters, results in a reliable estimation of 
the critical percolation exponents for the geometrical clusters in the limit of an infinite 
lattice size. The percolation exponents thus obtained, are also found to be consistent 
with those reported for the regular Ising model. These consistencies are explained 
in terms of the Fisher renormalization relations, which express the thermodynamic 
critical exponents of systems with annealed bond dilution in terms of those of the 
regular model system. 
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1. Introduction 

It is well known that the geometrical spin clusters (i.e., the clusters composed of 
the neighboring spins of the same sign), undergo a percolation transition at the 
thermodynamic critical temperature Tc in the two-dimensional Ising model [1]. This 
type of behavior is believed to be generally valid for a variety of two-dimensional 
critical models that undergo a continuous phase transition such as the g-state Potts 
model (g < 4) [2], which is a g-state generalization of the Ising model (g = 2) [3]. The 
formation of spanning spin clusters and their characteristic percolation exponents, can 
also be used to characterize the universality class of the corresponding thermal phase 
transition. Furthermore, at the critical point, the spanning cluster is a scale- invariant 
fractal object whose fractal dimensions uniquely specify the universality class of the 
associated continuous phase transition 

The critical behavior of a great number of statistical models in two spatial 
dimensions, has been investigated by the conformal field theory [5]. The conformal 
invariance property, refers to the invariance under coordinate transformations through 
which the angles between the crossing lines in the 2;-plane do not change. From 
this point of view, the cluster boundaries in the two-dimensional critical systems are 
considered as conformally invariant curves, and different characteristics such as their 
fractal dimensions are obtained. Indeed the fractal geometry is a useful mathematical 
tool for the characterization of a great many complex configurations. Self-similarity 
is the most important characteristic of such fractal objects. It is well known that the 
most widely studied statistical models in the condensed matter physics such as the 
Ising model, its g-state generalization the Potts model at criticality as well as the many 
critical geometrical phenomena exhibited by the various percolation models, consist of 
fractal lines H il [3 El El HO] . 

More recently, the spin cluster boundaries (interfaces) in the two-dimensional 
critical models, have been investigated rigorously using the method of Stochastic 
Loewner Evolution (SLE), in which the motion of a random walker along the cluster 
boundary in the upper-half complex plane in the continuum limit is specified by the 
Loewner dynamics 

dgtjz) _ 2 

dt gt{z) - ' ^ ^ 

which contains a Brownian term Q = v^-^* whose amplitude is given by the SLE 
parameter also known as the diffusivity [11]. The function gt maps a parametric 
curve 7t in the upper-half complex plane onto the real axis. Thus, given a real function 
C,t and using the initial condition g^ = z, the Loewner differential equation ([T]) determines 
the corresponding trace 7t in the upper-half ^-plane. The larger the diffusivity, k, the 
more is the deviation from a straight line. Indeed the nature of the SLE traces change 
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with the diffusivity: for the range of values < k < 4 — a range that includes k = 3 
characterizing the boundaries of the geometrical clusters (of like sign) in the critical 
two-dimensional Ising model — the SLE traces are nonintersecting simple curves, while 
for 4 < K < 8 the curves possess double points with possible self-touching (but no 
crossing), and for k > 8 they become space- filling [HI |T2l |T3]. Thus, the critical fractal 
dimension of the interfaces Dj introduced by this theory are model dependent. The 
relation between Dj and k, is given by = 1 + | [H] . The exact values are k = 3 and 
Dj = Y for the geometric spin cluster boundaries of the two-dimensional regular Ising 
model, as obtained analytically through various methods [15j. Like the thermodynamic 
critical exponents in statistical mechanics, k can divide different models into universality 
classes [13]. The main difference between k and the thermal critical exponents, is due 
to the method of definition — one being thermodynamic and the other geometrical in 
nature. Normally, the characterization of the universality class, requires a minimum 
of two thermodynamic critical exponents. We note that in the case of two-dimensional 
critical systems, a single parameter, k, appears to be sufficient to specify the universahty 
class. This may be attributed to the fact that of the two specifications of the model 
that significantly influence its critical behavior, i.e. space dimensions d and the order 
parameter dimensions n, one is held fixed at d = 2. Hence, the SLE diffusivity k alone 
can be used to specify the universality class of the two-dimensional critical systems. 

In the following we focus on the thermalized bond Ising model (TBIM) in two 
dimensions and investigate the behavior of its geometrical spin clusters (i.e., spin 
clusters of like sign) and their external boundaries (interfaces) at criticality. A Wolff 
single-cluster Monte Carlo algorithm is used to generate configurations at and near the 
criticality on square lattices, and a tie-breaking rule is used to identify non-intersecting 
geometrical cluster boundaries along the edges of the dual lattice. The rest of this 
paper is organized as follows. In section 2, the thermalized-bond Ising model is briefly 
reviewed, focussing on its thermodynamic critical behavior. The method of simulation 
and the finite size scaling procedures — employed to extrapolate the results obtained for 
finite lattices to the thermodynamic limit — are explained in section 3. The results are 
presented and discussed in section 4, and the paper is concluded with a summary in 
section 5. 

2. The model system 

The thermalized bond Ising model (TBIM), is a bond-diluted Ising model with a 
temperature dependent bond concentration, in which every covalent bond linking a 
nearest neighbor pair of atoms is allowed thermally induced electronic transitions 
between bonding and anti bonding electronic states [I6]. Hence, it can be regarded 
as containing annealed bond defects with a temperature dependent concentration. Each 
bond at every instant is characterized by a coupling constant Jij = 0, Jo, such that 
zero corresponds to a broken bond (anti-bonding electronic state), while Jq means an 
attractive coupling between the two atoms (bonding electronic state), as illustrated 
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schematically in figure [TJ Denoting the thermally averaged bond concentration by pb, 
(1 — Pb) must therefore represent the concentration of the broken bonds due to thermal 
excitations. To keep the analysis simple, the covalent bonds are treated as independent 
two-level systems with energy gap Jq, as sketched in figure [1] [16]. The ratio of the 
bonds to the broken bonds in equilibrium, is given by the ratio of the corresponding 
Boltzmann factors — Pb) = exp(/9Jo), or 

p, = l/(l + e-^^«) (2) 

where /3 = l/ksT is the reciprocal temperature and is the Boltzmann constant. The 
bond distribution function for the thermalized-bond model system is of the form 

^j.,(/3)=P6 5j,„Jo + (l-Pfe)5j.„o (3) 

where, pb is given by equation ([2]), and 5 denotes the Kronecker delta [17]. Hence, the 
Hamiltonian of the system can be formally defined by 

H = Y, J^S.S, (4) 

where Si = ±1 is an Ising spin, and the sum is over all nearest-neighbor pairs. 

It is well known that mapping the regular Ising model onto the equivalent correlated 
percolation problem introduces the bond probability 

P^ = l- e-^PJ (5) 

for the critical droplets (the Fortuin-Kastelyn clusters) of the regular Ising model [H HH] . 
The combination of equation ([2]) and equation ([5]), together with the choice Jq = 2 J, 
which is also the energy gap between parallel and antiparallel spins in the regular Ising 
model, results in a compound bond probability Ptbim = PrPb for our thermalized-bond 
model system: 

Ptbim = (1 - + e-^^') = tanh(/3J). (6) 

Equation ([6]) represents the bond probability used in our single-cluster update MC 
simulations. 

The thermodynamic critical behavior of the TBIM, has been studied before in two 
|19j . and three dimensions [16]. In two dimensions, the critical temperature is estimated 
to be Tc = 1.4897(3) [19], which is lower than the critical temperature of the regular Ising 
model. The lowering of the transition point is expected in the light of the annealed bond 
disorder present. However, the thermal critical exponents are found to be unchanged, 
within statistical errors. As for the three-dimensional TBIM, the thermal critical 
exponents are found to change consistent with the Fisher renormalization relations, 
as the specific heat exponent of the regular Ising model in three dimensions, — 0.11, 
is finite and positive. 
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Figure 1. Schematic illustration of the electronic energy states of the covalent bonds 
linking a chain of atoms. The energy gap between the bonding and the antibonding 
level is denoted by Jq. 



3. Simulation method and finite-size scaling 

As pointed out in the introduction, the critical behavior of the geometrical spin clusters 
and interfaces in the two-dimensional TBIM, is the main goal of this paper. For 
consistency with the postulates of SLE at T^, we have considered the model on strips 
of size X Ly, where the length of the strip is taken to be much larger than its 
width Ly = L with an aspect ratio L^/Ly = 8. The boundary conditions used for 
simulations, are fixed for the lower boundary (real axis), antiperiodic for the two sides, 
and free for the upper boundary of the system, as shown schematically in figure |2j 
Using a single-cluster update algorithm (Wolff's algorithm) [20] for the two-dimensional 
TBIM on square lattice, we generated equilibrium spin configurations at and near the 
critical point T^. A typical run consisted of several weeks of the CPU time on a single 
processor computer. Initially, the system was allowed 2x 10^ L equilibration Monte Carlo 
steps (MCS), and then the data points were accumulated by averaging over 2 x 10^ L 
configurations that contained a spanning cluster extending along the width of the strip 
L. Thus, L sets the appropriate length scale for the systems used, and the critical 
interfaces can be studied by the theory of SLE in the scaling limit. A turn right (or, 
alternatively, left) tie-breaking rule [21], is used on the square lattice as a procedure 
to identify the external perimeters (hulls) of the geometrical spin clusters without any 
self-intersection. We note that in this case the hulls and the external perimeters are the 
same. To identify the interfaces in the upper half complex plane, a walker moves along 
the edges of the dual lattice, starting from the origin as sketched in figure O At the first 
step of the walk, a spin (+) lies to the right of the walker (this direction is chosen to 
be the preferable direction). After arriving at each site on the dual lattice there are 3 
possibilities for the walker: it can cross any of the 3 nearest bonds of the original lattice. 
At the first step of selection, it chooses the bonds containing two different spins where 
crossing each of them leaves the spin (+) to the right and (— ) to the left of walker. The 
direction right or left are defined locally according to the orientation of the walker. If 
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there are still two possibilities for crossing, the walker chooses the bond which accords 
with the turn right tie-breaking rule. It turns toward the bond which is on its right-bond 
side with respect to its last direction in the last walk, if there is no selected bond to 
its right, it prefers to move straight on and if there is also no one there, it turns to its 
left. The procedure is repeated iteratively until the walker touches the upper boundary. 
The resulting interface is a unique one which has no self-intersection and never gets 
trapped [21]. Note that we just take the samples including a vertical spanning cluster 
in the y-direction. The fractal dimension of the interfaces at criticality, Dj, is obtained 
using the standard finite-size scaling procedure. The length of interfaces is related to 
the sample size as 

/ oc L^^ (7) 

Indeed the fractal dimension of the conformally invariant curves is provided by the SLE 
theory as 

= 1 + ^ (8) 

in which diffusivity k, as mentioned in the introduction, characterizes different 
universality classes, and so does Dj. For the regular Ising model, the diffusivity is 
believed to be k = 3, and thus Dj = ^ = 1.375. In addition, the fractal dimension of 
the spanning spin cluster at Tc, obeys the relation 

M oc L^^ (9) 

where M is the mass of the cluster, and obtained by counting all the nearest-neighbor 
positive (negative) spins to the right (left) of the SLE trace, as shown in figure [2j The 
exact value of Dc for the regular Ising model is Dc = ^ = 1.9479... [22] • Besides these, 
we find the winding angle variance through the winding angle function w{e) as defined 
by Wilson and Wieland [23]. For each edge on the dual lattice, there is a value for the 
winding function w{e) at that edge such that the winding angle at the neighboring edge 
e' is defined by w{e') = w{e) + the turning angle from e to e' measured in radians. 
It is shown that the variance of the winding angles, grows with the sample size like 

{e') = a+^\n{L). (10) 

2\ 



Thus, by plotting {6 ) versus ln(L), the slope gives a direct measure of k 
The finite-size scaling of the spanning cluster size, is of the form [21] 

M{L) = L^^M{L/i) (11) 

where the correlation length, also known as the connectedness length, behaves like 
,^ ~ (T — Tc)'"^ , and the scaling function M{x) tends to a constant as its argument 
goes to zero at Tc. Thus, the correlation length exponent of the geometrical clusters, 
is estimated by a value that results in a data collapse in a scaling plot L~^'=M(L) against 
L^/'^'^(/3//3c — 1). Among other percolation quantities of interest is the percolation 
strength P^{L) = M{L)/L'^, which is the probability that a site chosen at random 
belongs to the spanning cluster [21] . Poo plays the role of an order parameter for the 
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Figure 2. An schematic illustration of defining the domain boundaries in two- 
dimensional TBIM, on a square lattice. Figure shows the dual of the original square 
lattice including a spin configuration, with fixed boundary conditions for the bottom 
end, antiperiodic on sides, and free for the top end. The non-intersecting interface 
(shown by arrows) is generated using a turn right tie-breaking rule. 

percolation transition, and vanishes at the percolation threshold pc = tanh /3c, at a rate 
specified by an exponent fie defined by P^o oc {p~Vc)^'^ ■ The finite-size scaling relation 
for the strength of percolation, is as follows [2^ : 

P^{V) = L-'^«/^«P(L/0. (12) 

Thus, at the critical point T^, a log-log plot of Poo against L must be a straight line 
with a slope equal to the ratio —(3g/i^g- In the next section, we present our results for 
the two-dimensional TBIM. 

4. Results and discussion 

In this section we present and discuss our main results for the two-dimensional TBIM 
as obtained from simulations based on the methods pointed out in the previous section. 
We performed simulations for eight different system sizes L = 100, 150, 200, 250, 300, 
350, 400, and 450. Only the spin configurations including a vertical spanning cluster are 
considered for analysis, and the statistical errors were estimated by means of binning the 
accumulated data. As it appears in figure [3l the slope of a log-log plot of the spanning 
length / versus the system size L, results in a fractal dimension Dj = 1.373(8), which is 
equivalent to a k = 2.984(17) as given by equation (|8]). 

To confirm our results, we also measured the winding angle variance along the 
spanning contour by performing simulations for 10 different system sizes L = 30, 50, 
100, 150, 200, 250, 300, 350, 400, and 450. A plot of (6^) against L is shown in figure H 
A curve of the form ai + 02 ln(L) with parameters ai = —1.330(10) and 02 = 0.751(2) 
(where 02 = f), was least-squares fitted to these data. Furthermore, plotting (6*^) versus 
ln(L), results in a SLE parameter k = 3.004(9) as illustrated in the inset of figure HI 

A log-log plot of the spanning cluster mass M versus the system size L is shown in 
figure [5l The slope gives a fractal dimension Dc = 1.948(3). 

Using the finite-size scaling ansatz as given in equation (|TT|) . the correlation length 
exponent for the emerging spanning cluster, is estimated from a scaling plot of 
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Figure 3. A log-log plot of / against L for the two-dimensional TBIM. The slope of 
the straight line gives the fractal dimension of the external perimeters Dj — 1.373(8). 
The error bars are smaller or comparable with the symbol size. 
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100 200 300 400 

L 

Figure 4. A plot of {9"^) versus L for the two-dimensional TBIM. In the inset, the 
variance is in semi-logarithmic coordinates. The error bars are smaller or comparable 
with the symbol size. 

L~^''M{L) against V-/'^°{I3 / (3c — 1) as shown in figure El By varying i/g, and evaluating 
the quahty of the data collapse, our best estimate of the correlation length exponent 
for the geometrical clusters is vg = 1-01(2). Finally, the slope of the log-log plot of 
Poo versus L results in /3g = 0.051(3), as shown in figure [TJ which fits well into the 
hyperscaling relation for the percolation exponents in d spatial dimensions Dc = d— —. 

The obtained results for the two-dimensional TBIM are listed in table [T] for 
comparison with the analytical results of the regular Ising model. At criticality, the 
fractal dimension of the spin clusters Dc and interfaces Dj have been found to be 
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Figure 5. A log-log plot of M against L. The slope of the straight line gives the 
fractal dimension of the geometrical clusters at Tc, Dc — 1.948(3). 
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Figure 6. A scaling plot for the spanning cluster mass. By varying vg, and evaluating 
the quality of the data collapse, our best estimate of the correlation length exponent 
is VG = 1.01(2) 

Table 1. The critical exponents of the geometrical clusters in 2c? TBIM are compared 
with those obtained for 2d regular Ising model. 

VG Pg K £]_ £c 

TBIM 1.01(2) 0.051(3) 3.004(9) 1.373(8) 1.948(3) 

Ising Model 1.00 5/96 3 1.375 187/96 



consistent with the analytical results obtained for the regular Ising model, despite the 
temperature-dependent annealed bond dilution. 

It must be noted that although the geometrical clusters (i.e., the neighboring sites 
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Figure 7. A log-log plot of Poo against L. The slope of the straight line gives 
= 0.051(3) 

of the same spin sign) uniquely characterize the universahty class of the critical system, 
they do not, however, represent the critical droplets. Hence, the percolation critical 
exponents of the geometrical clusters, do not in general coincide with those of the 
corresponding thermal quantities. The critical droplets, are more precisely specified 
by the so-called Fortuin-Kastelyn (FK) clusters whose diffusivity parameter k has a 
duality relation with that of the geometric spin clusters [H [TBI EH] . The FK clusters can 
be obtained from the geometrical clusters through a random decimation of bonds by 
a suitable probability, in this case 1 — tanh(/3), and are therefore less compact. Thus, 
despite the correlation length exponent i^g = 1, we note that the geometrical clusters 
are too compact to represent the critical droplets, and the exponent [3g = 5/96 differs 
appreciably from the corresponding thermal exponent /3 = 1/8, associated with the 
vanishing of the magnetization order parameter at T^.. The value of the correlation length 
(also known as the connectivity length) exponent z/q = 1, for the geometrical clusters 
of the two-dimensional TBIM, is in excellent agreement with the values obtained from 
a real-space renormalization group analysis [H [7], high-temperature series expansion 
studies [26], and precision numerical simulations of the geometrical clusters of the 
standard two-dimensional Ising model [21 127]. This result, however, opens a question 
about a near perfect collapse onto a universal function for the same data obtained for 
the regular Ising model, but with a different exponent ^ [2T] . 

As can be seen from table [1], within the statistical uncertainty, the value of the SLE 
parameter k, the fractal dimensions, and the percolation exponents of the geometrical 
spin clusters of the two-dimensional TBIM, are in excellent agreement with those of 
the corresponding regular Ising model. These results agree well with an earlier study 
of the thermodynamic critical behavior of the two-dimensional TBIM, which places the 
model system in the universality class of the standard two-dimensional Ising model 
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|19] . and the Fisher renormahzation relations, which assert that annealed bond dilution 
can only change the critical exponents if the specific heat exponent of the regular 
model is positive (a^ > 0) [28]. The two-dimensional regular Ising model, however, is 
characterized by a logarithmic divergence of the specific heat, = 0, and the exponents 
remain unchanged. 

As for the fractal behavior of the geometrical clusters away from the criticality, we 
note that for all temperatures below the transition point T < Tc (or p > pc), must 
equal the space dimensions d = 2, otherwise the percolation strength would vanish in 
the thermodynamic limit of an infinite lattice. At {p = Pc), the geometrical clusters 
of the two-dimensional TBIM, are fractal with a fractal dimension Dc = 1.948(3) (< 2), 
thus rendering the percolation strength zero at the critical point, as expected. Hence, 
one expects the fractal dimension of the geometrical clusters to change discontinuously 
from Dc = d = 2 for T < Tc, to Dc = 1.948(3) at Tc. This expectation is validated 
by the data of reference [21], where the so-called 'effective' fractal dimensions undergo 
a sharp crossover at Tc. We believe that the crossover is a finite-size effect, and a 
remnant of the discontinuity at Tc in the thermodynamic limit. As for the temperatures 
above the critical point T > Tc, there are no spanning geometrical clusters in the 
thermodynamic limit and the procedures used here become inapplicable. However, other 
standard procedures such as the box counting method may be employed to investigate 
the fractal behavior of the geometrical clusters above Tc and within a region of linear 
size of the order of the finite correlation length ^. 

5. Summary 

The fractal behavior of the geometrical spin clusters, are obtained for the thermalized 
bond Ising model in two dimensions. For this purpose, a modified Wolff single-cluster 
Monte Carlo simulation is used to generate equilibrium spin configurations on square 
lattices in the critical region. The obtained values for the fractal dimensions of the 
spanning geometrical clusters, Dc, and that of their interfaces, Dj, are in perfect 
agreement with those reported for the regular Ising model. The variance of the 
winding angles results in a value k = 3.004(9) for the SLE parameter, thus placing 
it in the universality class of the regular Ising model. Furthermore, the percolation 
exponents of the geometrical spin clusters at Tc, are found to be consistent with those 
reported for the regular Ising model. These consistencies are explained in terms of the 
Fisher renormahzation relations, which express the thermodynamic critical exponents 
of systems with annealed bond dilution in terms of those of the regular model system. 
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